!	** eamfit **
!	Fit (force-match) EAM potential parameters to produce reference forces

	module rijdata
	implicit none
	! Convenience module containing allocatable rij() and eamtypes arrays, and natoms
	! Why? The cost function is a generic wrapper called by the Simplex minimiser, and needs to
	! be able to call eam_forces() with the necessary data.
	integer :: natoms = 0
	real*8, allocatable :: rij(:,:,:)
	integer, allocatable :: eamtypes(:)
	! Variables for Simplex Simulated Annealing fit
	integer :: fit_nframes = 0, frame
	real*8 :: fit_sose
	real*8, allocatable :: fit_r(:,:,:), fit_f(:,:,:), calc_f(:,:)
	real*8, allocatable :: fit_cell(:,:), fit_icell(:,:), fit_opteam(:,:,:)

	contains

	subroutine allocate_rijdata
	implicit none
	! Test for existing alloc, and deallocate if necessary
	if (allocated(eamtypes)) deallocate(eamtypes, rij, fit_r, fit_f, fit_cell, fit_icell, calc_f)
	allocate(eamtypes(natoms), rij(natoms,natoms,4), fit_r(fit_nframes,natoms,3), fit_f(fit_nframes,natoms,3), fit_cell(fit_nframes,9), fit_icell(fit_nframes,9), calc_f(natoms,3))
	end subroutine allocate_rijdata

	end module rijdata

	program eamfit
	use parse; use eamdata; use simplex; use cellutils; use mdroutines; use random; use rijdata
	implicit none
	! General variables
	character*80 :: infile
	character*20 :: temp
	integer :: n, m, o, i, iargc, nargs, seed = 9281750
	! FCC cell definition
	integer :: fcc_nx, fcc_ny, fcc_nz
	real*8 :: fcc_alpha, fcc_jiggle = -1.0d0
	logical :: fcc_defined = .false.
	! Variables for Molecular Dynamics simulation
	integer :: step, md_writefreq = 1
	logical :: md_usevelocityverlet = .false.
	! Variables for numerical force calculation test
	character*12 :: test_errorstring
	real*8 :: test_numforces_phi(3), test_numforces_F(3), test_tempr(3), test_delta = 0.01d0, test_energy_phi(2), test_energy_F(2)
	real*8, allocatable :: test_forces_phi(:,:), test_forces_F(:,:)
	! Variables for energy-volume curve calculation
	real*8 :: ev_delta
	integer :: ev_nsteps
	! Program flags
	logical :: perform_fit = .false., perform_md = .false., perform_forcetest = .false., perform_energytest = .false.
	logical :: perform_ev = .false.
	! External function declarations
	real*8, external :: cost, energy_and_forces

	nargs = iargc()
	if (nargs.eq.0) then
	  write(0,"(a)") "Usage : eamfit <datafile> [options...]"
	  write(0,"(a)") "           # SSA fit control"
	  write(0,"(a)") "               [-fit npasses nloops ncycles] [-fix eamtype parameterid]"
	  write(0,"(a)") "           # Molecular dynamics control"
	  write(0,"(a)") "               [-md nsteps timestep temperature]"
	  write(0,"(a)") "               [-writefreq steps] [-vv]"
	  write(0,"(a)") "           # FCC lattice generation"
	  write(0,"(a)") "               [-fcc alpha nx ny nz] [-jiggle delta]"
	  write(0,"(a)") "           # Simple energy/force (test) options"
	  write(0,"(a)") "               [-force deltar] [-energy] [-ev delta nsteps]"
	  write(0,"(a)") "           # General options"
	  write(0,"(a)") "               [-seed i] [-nocut]"
	  stop
	end if
	call getarg(1,infile)
	
	! Open and read input file
	write(6,*) "Reading input file..."
	open(unit=11,file=infile,form='formatted',status='old',err=20)
	! First, number of atom types / EAM potentials defined, and their functional form
	if (.not.readline(11)) goto 20
	if (argi(1).gt.1) stop "Only one EAM type is allowed for now."
	eam_form = argi(2)
	call alloc_eam(argi(1))
	eamnames = ""
	eam = 0.0
	eamflags = 1
	if (eam_form.eq.1) write(6,"(24x,a8,1x,12a8)") "Name    ", "  Mass  ", "   r0   ", "   e0   ", "   A    ", " alpha  ", "  beta  ", "  rho0  ", "   Zi   "
	if (eam_form.eq.2) write(6,"(24x,a8,1x,14a8)") "Name    ", "  Mass  ", "   Ec   ", "  phi0  ", "   r0   ", " alpha  ", "  beta  ", " gamma  ", " delta  ", "   c0   ", "   c1   ","   c2   ", "   c3   ", "   c4   ", "   c5   "
	do n=1,neam
	  if (.not.readline(11)) goto 20
	  ! Format of line is : "name  mass(g)  r0(Ang)  e0(eV)   A   alpha(Ang)  beta  rho0  Zi"
	  eamnames(n) = arg(1)(1:8)
	  eammass(n) = argr(2)
	  do m=1,NPARAMS(eam_form)
	    eam(n,m) = argr(m+2)
	  end do
	  write(6,"(a20,i2,2xa8,14f8.4)") "  -- Read EAM data :", n, eamnames(n), eammass(n), eam(n,1:NPARAMS(eam_form))
	end do
	inputeam = eam
	! Next, number of atoms in each frame, and then number of frames to expect in file
	if (.not.readline(11)) then
	  write(6,"(a,i4)") "Error reading number of atoms/framesfrom file"
	  goto 20
	end if
	natoms = argi(1)
	fit_nframes = argi(2)
	if ((natoms.gt.0).and.(fit_nframes.gt.0)) call allocate_rijdata()
	eamtypes = 0
	write(6,"(a,i4,a,i3,a)") "Expecting ", fit_nframes, " frames containing ", natoms, " atoms each"
	do frame=1,fit_nframes
	  ! First data is unit cell specification as nine real numbers on one line
	  if (.not.readline(11)) then
	    write(6,"(a,i4)") "Error reading cell for frame ",frame
	    goto 20
	  end if
	  do m=1,9
	    fit_cell(frame,m) = argr(m)
	  end do
	  fit_icell(frame,:) = fit_cell(frame,:)
	  call calc_icell(fit_icell(frame,:))
	  ! Then follow 'natoms' lines of the format :  "eamname  x  y  x  fx  fy  fz"
	  do i=1,natoms
	    if (.not.readline(11)) then
	      write(6,"(a,i4,a,i4)") "Error reading atom ",i," for frame ",frame
	      goto 20
	    end if
	    ! Determine eam type to use for this atom (first frame only)
	    if (frame.eq.1) then
	      do m=1,neam
		if (arg(1).eq.eamnames(m)) eamtypes(i) = m
	      end do
	      if (eamtypes(i).eq.0) then
		write(6,"(a,i,a,a)") "Failed to find corresponding type for atom ",i, "name = ",arg(1)
		stop
	      end if
	    end if
	    fit_r(frame,i,1) = argr(2)
	    fit_r(frame,i,2) = argr(3)
	    fit_r(frame,i,3) = argr(4)
	    fit_f(frame,i,1) = argr(5)
	    fit_f(frame,i,2) = argr(6)
	    fit_f(frame,i,3) = argr(7)
	  end do
	  write(6,"(a,i4,a,i4)") "Successfully read frame ",frame," of ",fit_nframes
	end do
	goto 25
20	stop "Error reading input file."
25	close(11)

	! Now handle command-line arguments
	n = 1
	do
	  n = n + 1; if (n.gt.nargs) exit
	  call getarg(n,temp)
	  select case (temp)
	    ! SSA Switches
	    case ("-fit")
	      n = n + 1; call getarg(n,temp); read(temp,"(I6)") simplex_npasses
	      n = n + 1; call getarg(n,temp); read(temp,"(I6)") simplex_nloops
	      n = n + 1; call getarg(n,temp); read(temp,"(I6)") simplex_ncycles
	      perform_fit = .true.
	    case ("-fix")
	      n = n + 1; call getarg(n,temp); read(temp,"(I6)") i
	      n = n + 1; call getarg(n,temp); read(temp,"(I6)") m
	      if ((i.lt.1).or.(i.gt.neam)) stop "Argument error: EAM type index is out of range for the number of defined EAM types."
	      if (m.gt.NPARAMS(eam_form)) stop "Argument error: EAM parameter index is out of range."
	      eamflags(i,m) = 0
	    ! MD Switches
	    case ("-md")
	      n = n + 1; call getarg(n,temp); read(temp,"(I6)") md_nsteps
	      n = n + 1; call getarg(n,temp); read(temp,"(f12.6)") md_timestep
	      n = n + 1; call getarg(n,temp); read(temp,"(f12.6)") md_Ttarget
	      perform_md = .true.
	    case ("-writefreq")
	      n = n + 1; call getarg(n,temp); read(temp,"(i9)") md_writefreq
	    case ("-vv")
	      md_usevelocityverlet = .true.
	    ! FCC Cell Specification
	    case ("-fcc")
	      n = n + 1; call getarg(n,temp); read(temp,"(f12.6)") fcc_alpha
	      n = n + 1; call getarg(n,temp); read(temp,"(i6)") fcc_nx
	      n = n + 1; call getarg(n,temp); read(temp,"(i6)") fcc_ny
	      n = n + 1; call getarg(n,temp); read(temp,"(i6)") fcc_nz
	      fcc_defined = .true.
	    case ("-jiggle")
	      n = n + 1; call getarg(n,temp); read(temp,"(f12.6)") fcc_jiggle
	    ! Energy/Force Calculation Switches
	    case ("-force")
	      n = n + 1; call getarg(n,temp); read(temp,"(f12.6)") test_delta
	      perform_forcetest = .true.
	      perform_energytest = .true.
	    case ("-energy")
	      perform_energytest = .true.
	    case ("-ev")
	      n = n + 1; call getarg(n,temp); read(temp,"(f12.6)") ev_delta
	      n = n + 1; call getarg(n,temp); read(temp,"(i6)") ev_nsteps
	      perform_ev = .true.
	    ! General Switches
	    case ("-seed")
	      n = n + 1; call getarg(n,temp); read(temp,"(I6)") seed
	    case ("-nocut")
	      eam_cutoffs = .false.
	    case default
	      write(0,"(a,a)") "Unrecognised command line option:",temp
	      stop
	  end select
	end do

	! Write info!
	write(6,"(A,I9)") "Random seed is ",seed
	call initrandom(seed)

	if (perform_fit) then
	  write(6,"(/,a)") "############################"
	  write(6,"(a)") "# Beginning Simplex SA Run #"
	  write(6,"(a,/)") "############################"

	  write(6,"(a,i)") "Total number of parameters to adjust is", sum(eamflags)
	  ! Construct array of initial alpha values (fit parameters) - all parameters are added, regardless of fit flags (this should be handled in the cost function)
	  ! Allocate alpha array and set nalpha - note, this includes *all* parameters, regardless of the corresponding fit flag
	  ! Fixing of parameters should be handled in the 'cost()' function when translating the passed alpha array back into usable potential parameters
	  simplex_nalpha = neam*NPARAMS(eam_form)
	  call simplex_createarrays()
	  ! Copy starting parameters to alpha array
	  simplex_alpha(:) = reshape( eam, (/ neam*4 /) )

	  ! Loop over and fit to frames
	  fit_opteam = 0.0
	  do frame=1,fit_nframes
	    write(6,"(a,i4)") " .. Moving to frame ", frame
	    call calc_rij(natoms, fit_cell(frame,1:9), fit_icell(frame,1:9), fit_r(frame,1:natoms,1:3), rij)
	    write(6,*) " .. Initial cost = ", cost(simplex_nalpha, simplex_alpha)
	    ! Perform minimisation (if there is at least one variable parameter)
	    if (sum(eamflags).ne.0) then
	      write(6,"(a,i,a,i,a,i,a)") "Beginning ", simplex_npasses, " minimisations comprising ", simplex_nloops, " loops of ", simplex_ncycles, " SimplexSA cycles..."
	      do n=1,simplex_npasses
	        fit_sose = simplex_minimise()
	        fit_opteam(frame,:,:) = reshape( simplex_alpha, (/ neam,NPARAMS(eam_form) /) )
	      end do
	    end if
	    write(6,*) " .. Final Cost = ", fit_sose
	  end do

	  write(6,*) "Optimized parameters:"
! 	  write(6,*) "       X           Y           Z           Q"	
	  write(6,*) eam
	  !do f=1,fit_nframes
	  !  write(6,"(a,4f12.6)")  points(n,1:3), charges(n)
	  !end do
	end if

	! All other run types might (or need to) use a defined FCC lattice...
	if (fcc_defined) then
	  call md_createfcc(fcc_alpha,fcc_nx,fcc_ny,fcc_nz,eammass(1),eammass(1),eammass(1),eammass(1),.false.)
	  ! Lastly, reallocate working arrays to reflect new size of system
	  natoms = md_natoms
	  fit_nframes = 1
	  call allocate_rijdata()
	  eamtypes = 1
	  md_atomnames = eamnames(eamtypes(1))
	  ! Jiggle initial positions?
	  if (fcc_jiggle.gt.0.0d0) then
	    do n=1,natoms
	      test_tempr = (/ ran()-0.5d0, ran()-0.5d0, ran()-0.5d0 /)
	      test_tempr = (test_tempr / sqrt(sum(test_tempr*test_tempr))) * fcc_jiggle
	      md_r(n,:) = md_r(n,:) + test_tempr(:)
	    end do
	  end if
	  ! Copy data to fit_* arrays - all but the MD code require it to be there
	  fit_r(1,:,:) = md_r(:,:)
	  fit_cell(1,:) = md_cell(:)
	  fit_icell(1,:) = md_icell(:)
	end if

	if (perform_ev) then
	  write(6,"(/,a)") "###############################"
	  write(6,"(a)") "# Beginning Energy/Volume Run #"
	  write(6,"(a,/)") "###############################"
	  ! Starting configuration a defined crystal lattice
	  if (.not.fcc_defined) stop "An initial FCC lattice must be provided with '-fcc'"
	  ! Print EV run parameters
	  write(6,"(a50,es10.4)") "Initial lattice parameter alpha (Angstroms) = ", fcc_alpha
	  write(6,"(a50,i6)") "Number of steps to take = ", ev_nsteps
	  write(6,"(a50,es10.4)") "Stepsize for lattice parameter (Angstroms) = ", ev_delta
	  write(6,"(a50,es10.4)") "Final lattice parameter alpha (Angstroms) = ", fcc_alpha + ev_nsteps*ev_delta
	  ! Calculate EV curve
	  open(unit=15,file="ev.out",form='formatted',status='replace')
	  write(6,"(a4,a8,4(a12,2x))") "Step", "Alpha", "Volume", "V/Atom", "Etot(eV)", "Epair(eV)", "Eemb(eV)"
	  write(15,"('#',a4,a8,5(a12,2x))") "Step", "Alpha", "Volume", "V/Atom", "Etot(eV)", "Epair(eV)", "Eemb(eV)"
	  do n=0,ev_nsteps
	    call calc_rij(natoms, md_cell, md_icell, md_r, rij)
	    call eam_energy_and_forces(natoms, eamtypes, rij)
	    write(6,"(i4,f8.4,5(es12.4,2x))") n, fcc_alpha, volume(md_cell), volume(md_cell)/md_natoms, eam_energy_total, eam_energy_phi, eam_energy_F
	    write(15,"(i4,f8.4,5(es12.4,2x))") n, fcc_alpha, volume(md_cell), volume(md_cell)/md_natoms, eam_energy_total, eam_energy_phi, eam_energy_F
	    fcc_alpha = fcc_alpha + ev_delta
	    call md_createfcc(fcc_alpha,fcc_nx,fcc_ny,fcc_nz,eammass(1),eammass(1),eammass(1),eammass(1),.true.)
	  end do
	  close(15)
	end if

	if (perform_md) then
	  write(6,"(/,a)") "####################"
	  write(6,"(a)") "# Beginning MD Run #"
	  write(6,"(a,/)") "####################"
	  ! Starting configuration is either first frame of trajectory supplied, or a defined crystal lattice (created earlier)
	  if (.not.fcc_defined) then
	    ! Allocate necessary variables
	    call md_createempty(natoms)
	    ! Assign atomic masses and atom names
	    do i=1,md_natoms
	      md_mass(i) = eammass(eamtypes(i))
	      md_atomnames(i) = eamnames(eamtypes(1))
	    end do
	    ! Grab initial particle positions and cell from first frame in supplied trajectory
	    md_r(:,:) = fit_r(1,:,:)
	    md_cell = fit_cell(1,:)
	    md_icell = md_cell
	    call calc_icell(md_icell)
	  end if
	  ! Print MD run parameters
	  write(6,"(a50,i6)") "Number of steps = ", md_nsteps
	  write(6,"(a50,es10.4)") "Simulation timestep (ps) = ", md_timestep
	  write(6,"(a50,es10.4)") "Initial temperature (Kelvin) = ", md_Ttarget	
	  write(6,"(a50,i6)") "Trajectory/energy output frequency = ", md_writefreq
	  if (md_usevelocityverlet) write(6,"(a50)") "Using velocity Verlet integrator"

	  ! Open output file
	  open(unit=15,file="md.out",form='formatted',status='replace')
	  ! Set initial random velocities
	  call md_initialise()
	  ! Write starting configuration
	  call md_writeframe(.true.)
	  ! Calculate potential energy and forces, and then kinetic energy and temperature
	  call calc_rij(natoms, md_cell, md_icell, md_r, rij)
	  call eam_energy_and_forces(natoms, eamtypes, rij)
	  md_pe = eam_energy_total
	  call md_calc_kinetic()
	  write(6,"(a6,6(a12,2x))") "", "", "", "Energies, eV   "
	  write(6,"(a6,8(a12,2x))") "Step", "T, K    ", "Total   ","Kinetic   ", "Potential", "P.E. / Atom", "P.E. (emb)", "P.E. (pair)", "COM Velocity"
	  write(6,"(i5,1x,7(es12.5,2x),3f6.3)") 0, md_Tcalc, md_ke+md_pe, md_ke, md_pe, md_pe/natoms, eam_energy_F, eam_energy_phi, md_sumv
	  write(15,"('#',a6,8(a12,2x))") "Step", "T, K    ", "Total   ","Kinetic   ", "Potential", "P.E. / Atom", "P.E. (emb)", "P.E. (pair)", "COM Velocity"
	  write(15,"(i5,1x,7(es12.5,2x),3f6.3)") 0, md_Tcalc, md_ke+md_pe, md_ke, md_pe, md_pe/natoms, eam_energy_F, eam_energy_phi, md_sumv
	  do md_step=1,md_nsteps
	    call calc_rij(natoms, md_cell, md_icell, md_r, rij)
	    call eam_energy_and_forces(natoms, eamtypes, rij)
	    md_pe = eam_energy_total
	    if (md_usevelocityverlet) then
	      call md_velocity_verlet(eam_forces_total)
	    else
	      call md_verlet(eam_forces_total)
	    end if
	    call md_calc_kinetic()
	    if (mod(md_step,md_writefreq).eq.0) then
	      write(6,"(i5,1x,7(es12.5,2x),3f6.3)") md_step, md_Tcalc, md_ke+md_pe, md_ke, md_pe, md_pe/natoms, eam_energy_F, eam_energy_phi, md_sumv
	      write(15,"(i5,1x,7(es12.5,2x),3f6.3)") md_step, md_Tcalc, md_ke+md_pe, md_ke, md_pe, md_pe/natoms, eam_energy_F, eam_energy_phi, md_sumv
	      call md_writeframe(.false.)
	    end if
	  end do
	  close(15)
	  ! Free allocated arrays
	  call md_freearrays()
	end if

	if (perform_forcetest.or.perform_energytest) then
	  write(6,"(/,a)") "###############################"
	  write(6,"(a)") "# Beginning Energy/Force Test #"
	  write(6,"(a,/)") "###############################"
	  ! Allocate necessary variables
	  allocate(test_forces_phi(natoms,3), test_forces_F(natoms,3))
	  do frame = 1,fit_nframes
	    write(6,"(/,'Frame ',i4)") frame
	    call calc_rij(natoms, fit_cell(frame,1:9), fit_icell(frame,1:9), fit_r(frame,1:natoms,1:3), rij)
	    call eam_energy_and_forces(natoms, eamtypes, rij)
	    write(6,"(a13,e12.5,a,3x,e12.5,a,e12.5,a,e12.5)") "Total P.E. = ", eam_energy_total, " eV, ", eam_energy_phi, " eV (pair), ", eam_energy_F, " eV (emb), A (angs) = ", fit_cell(frame,1) 

	    if (.not.perform_forcetest) cycle
	    ! Copy current, calculated forces into temporary arrays
	    test_forces_phi = eam_forces_phi
	    test_forces_F = eam_forces_F
	    write(6,"(/,12x,(7x,a12,8x),2(x,a25,x),(7x,a12,8x))") "Angstroms", "eV/Angstrom (Analytic)", "eV/Angstrom (Numeric)", "     Delta      "
	    write(6,"(a9,4x,6(3x,a2,4x),3(3x,a3,3x),3(x,a6,3x))") "ID / Atom", "rx", "ry", "rz", "fx", "fy", "fz", "nfx", "nfy", "nfz", "fx-nfx", "fy-nfy", "fz-nfz"
	    do i=1,natoms
	      test_tempr(1:3) = fit_r(frame,i,1:3)
	      ! Calculate numerical forces for each atom in each direction
	      do n=1,3
		fit_r(frame,i,n) = test_tempr(n) - test_delta
		call calc_rij(natoms, fit_cell(frame,1:9), fit_icell(frame,1:9), fit_r(frame,1:natoms,1:3), rij)
		call eam_energy_and_forces(natoms, eamtypes, rij)
		test_energy_phi(1) = eam_energy_phi
		test_energy_F(1) = eam_energy_F
		fit_r(frame,i,n) = test_tempr(n) + test_delta
		call calc_rij(natoms, fit_cell(frame,1:9), fit_icell(frame,1:9), fit_r(frame,1:natoms,1:3), rij)
		call eam_energy_and_forces(natoms, eamtypes, rij)
		test_energy_phi(2) = eam_energy_phi
		test_energy_F(2) = eam_energy_F
		fit_r(frame,i,n) = test_tempr(n)
		! Calculate numerical forces
		test_numforces_phi(n) = (test_energy_phi(1)-test_energy_phi(2)) / (2.0d0*test_delta)
		test_numforces_F(n) = (test_energy_F(1)-test_energy_F(2)) / (2.0d0*test_delta)
	      end do
	      ! Write force comparison
	      test_errorstring = ""
	      if (sum((test_forces_phi(i,1:3)-test_numforces_phi(1:3))**2).gt.1d-5) test_errorstring = "!!! Error"
	      write(6,"(i4,x,a4,2x,9f9.4,3es10.2,2x,a3,2x,a12)") i,eamnames(eamtypes(i)), fit_r(frame,i,1:3), test_forces_phi(i,1:3), test_numforces_phi(1:3), test_forces_phi(i,1:3)-test_numforces_phi(1:3), "phi", test_errorstring
	      if (sum((test_forces_phi(i,1:3)-test_numforces_phi(1:3))**2).gt.1d-5) write(6,"(92x,3es10.2,7x,a12)") test_forces_phi(i,1:3)/test_numforces_phi(1:3), "(ratios)"
	      test_errorstring = ""
	      if (sum((test_forces_F(i,1:3)-test_numforces_F(1:3))**2).gt.1d-5) test_errorstring = "!!! Error"
	      write(6,"(38x,6f9.4,3es10.2,2x,a3,2x,a12)") test_forces_F(i,1:3), test_numforces_F(1:3), test_forces_F(i,1:3)-test_numforces_F(1:3), "F", test_errorstring
	      if (sum((test_forces_F(i,1:3)-test_numforces_F(1:3))**2).gt.1d-5) write(6,"(92x,3es10.2,7x,a12)") test_forces_F(i,1:3)/test_numforces_F(1:3), "(ratios)"
	    end do
	  end do
	end if
	
	end program eamfit

	subroutine calc_rij(natoms, cell, icell, r, rij)
	use cellutils
	implicit none
	integer :: i,j,natoms
	real*8 :: cell(9), icell(9), r(natoms,3), rij(natoms,natoms,4), vecij(3)
	! Calculate the distance/vector matrix between atoms for speed and simplicity
	rij = 0.0
	do i=1,natoms-1
	  do j=i+1,natoms
	    call mimvec(cell, icell, r(j,1), r(j,2), r(j,3), r(i,1), r(i,2), r(i,3), vecij(1), vecij(2), vecij(3))
	    rij(i,j,1:3) = vecij(1:3)
	    rij(i,j,4) = sqrt(vecij(1)*vecij(1)+vecij(2)*vecij(2)+vecij(3)*vecij(3))
	    rij(j,i,1:3) = -vecij(1:3)
	    rij(j,i,4) = rij(i,j,4)
! 	    write(0,"(2i4,4f12.6)") i,j,rij(i,j,:)
	  end do
	end do
	end subroutine calc_rij

	real*8 function cost(nalpha,alphaparams)
	use rijdata; use eamdata
	implicit none
	integer :: nalpha, n, m, i, j
	real*8, intent(inout) :: alphaparams(nalpha)
	real*8 :: energy, calcf(natoms,3)

	! Copy current fit parameter data back into local data structure
	! Overwrite any parameters with their original input values if fitflag = 0
	eam(:,:) = reshape( alphaparams, (/ neam,NPARAMS(eam_form) /) )
	do n=1,neam
	  do m=1,NPARAMS(eam_form)
	    if (eamflags(n,m).eq.0) eam(n,m) = inputeam(n,m)
	  end do
	end do

	call eam_energy_and_forces(natoms, eamtypes, rij)

	! Calculate cost function
	cost = 0.0
	do n=1,natoms
	  cost = cost + (fit_f(frame,n,1)-eam_forces_total(n,1)) * (fit_f(frame,n,1)-eam_forces_total(n,1))
	  cost = cost + (fit_f(frame,n,2)-eam_forces_total(n,2)) * (fit_f(frame,n,2)-eam_forces_total(n,2))
	  cost = cost + (fit_f(frame,n,3)-eam_forces_total(n,3)) * (fit_f(frame,n,3)-eam_forces_total(n,3))
	end do
	! Copy (fixed) parameters back to alpha array
	alphaparams(:) = reshape( eam, (/ neam*4 /) )
	end function cost
